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Abstract 


A two-dimensional network model for porous media was developed to simulate the 
transport of bacteria through saturated porous media. The model was a regular square 
network consisting of pore throats or tubes with randomly distributed radii from Weibull 
distribution. The length of the pore throat was calculated by equating the porosity of the 
porous media with that of the network. The network thus obtained simulated the reported 
results of drainage capillary pressure curve for porous media. Given the heads at the inlet 
and outlet, head at each node was computed by repeated applications of Darcy-Weisbach 
equation in the network and solving the system of non-linear equations using the Newton- 
Raphson method. The hydraulic model was validated by hand-computed values in simple 
networks. The particle transport model was incorporated by repeated application of 
Iwasaki equation to each tube in the network with the assumption that the filter coefficient 
is constant in a small tube during a small time step. Using the mass balance equation, 
deposition in the tubes was also computed. The filter coefficients and velocities in the 
tubes were updated at each time step using the deposition values. The transport model was 
checked for consistency by simulating for conservative tracer. The simulation results for a 
non-conservative particle (bacteria) transport shows similar qualitative trends as reported 
for bacterial transport column experiments. The effects of various distribution parameters 
and grain sizes on the pore throat (tube) radii distribution are discussed. The model was 
then used to study the effects of various parameters such as, porosity, flow rate, adhesion 
efficiency (a), maximum specific deposit ( o max ) etc. on breakthrough curves. The adhesion 
efficiency (a) was found to be a critical parameter that can significantly alter the transport. 
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Chapter 1 Introduction 

Chapter 1 

Introduction 

Frequently ground water is used without treatment because of the perceived filtering 
action of solution flow through porous media. It is now recognized, however, that 
bacteria and viruses can travel considerable distances in aquifers and saturated soils, 
thus posing a contamination threat to surface waters and well waters (Bouwer, 1984). 
Historically, bacterial contamination of ground water can be associated with the age-old 
practices of traditional septic tank disposal systems and land application of sewage 
sludge (Gerba and Bitton, 1984). However the concern over microbial quality of 
ground water supplies is fairly recent. A number of field studies have indicated that 
many waterborne epidemics owe their origin to bacterial contamination of ground 
water (Wilson et al., 1983, Romero, 1970). It is well established that, microbial 
contamination occurs when human wastes enter the soil, since microbes can travel long 
distances through aquifers under appropriate conditions (Corapcioglu and Haridas, 
1984). A large number of outbreaks of water-borne disease have been attributed to 
contaminated groundwater (Gerba and Bitton, 1984). Control of microbial transport 
and retention in porous media is also required for effective and safe bioremediation 
(Bouwer, 1993). Thus the study of bacterial transport is of paramount importance, in 
order to determine whether the bacteria, introduced to the subsurface reach the point of 
interest. 

Bacterial mobility through a porous medium depends on a number of processes such as, 
advection, dispersion, attachment, detachment, etc (Smith et ah, 1985; Harvey et ah, 
1993). Body forces arising due to gravity, buoyancy, hydrodynamic drag, Van-der- 
Waals interaction, double layer interaction and hydrophobic interaction affect the 
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transport, attachment and detachment of bacteria (Rajagopalan and Tien, 1976). 
Natural heterogeneity of porous media adds further complexity to the transport process 
(Fontes et al., 1991; Harvey and Garabedian 1991). As a result, mechanistic models for 
complete bacterial transport that takes into account the effects of individual processes 
and forces are non-existent. Large-scale transport models often consider bacterial 
retention on the porous media either as an equilibrium process or as a first order kinetic 
process, devoid of any mechanistic description of individual forces (Corapcioglu and 
Haridas, 1984). Attempts to incorporate the mechanisms of bacterial retention on 
porous media typically consider the process to be a combination of transfer of bacteria 
from bulk water to the collector surface (mass transfer efficiency) and adhesion 
(collision efficiency). The mass transfer efficiency was simulated (Rajagopalan and 
Tien, 1976; Basu, 2001) using the sphcre-in-cell porous media model (Happel, 1958). 
This lead to initial collection efficiency when the collision efficiency was assumed to 
be unity. Validation of these models requires computation of initial collection 
efficiencies from the breakthrough curves obtained from laboratory column 
experiments for comparison with the model simulations. Computation of the initial 
collection efficiency from discrete time series (experimental breakthrough curve) can 
be cumbersome and introduce transformation errors (Rijnaarts et al. 1996). 

A model that incorporates the physical processes as well as capable of directly 
simulating the breakthrough curve will be helpful in predicting laboratory column 
experiments results under changing physical parameters. Such a model will be a robust 
tool for design of experiments and prediction of qualitative effect of different physical 
parameters on bacterial transport. Goal of the present work is to develop a 2-D network 
model for bacterial transport through porous media under saturated flow condition that 
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qualitatively simulates the changes in breakthrough curves obtained from laboratory 
column experiments under changing physical parameters. 

The thesis has been organized as follows: 

Chapter 2 gives the status of the literature on the related areas. 

Chapter 3 states the scope of the work. 

Chapter 4 outlines the development of the model. 

Chapter 5 describes the results obtained from the simulations and compares with the 
reported experimental results. 

Chapter 6 summarizes the work and states the key conclusions. 

Chapter 7 lists the scope for future work. 
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2.1 Modeling of Micro-structure and Fluid Flow through 
Porous Media 

The earliest studies of fluid flow through porous media used the sphere pack as model. 
Kozney (1927) and Carman (1938) derived an equation, which related permeability of 
porous media to porosity and internal surface area of sphere pack. Rapport and Leas 
(1952) reduced the porous medium to an equivalent bundle of tubes through many 
assumptions. Many frequency distributions were proposed that related the properties of 
porous media to the radii of the equivalent bundle of tubes (Childs and Collis-George, 
1948; Gates and Tempelaar Leitz, 1950; Fatt and Dykstra, 1951; Purcell, 1949; 
Burdine, et al., 1950). Several obvious characteristics of real porous media are not 
present in the bundle of tubes model. The major weakness of the bundle of tubes 
model lies in the absence of cross-connections between the tubes. An examination of 
thin sections of sandstone indicated that the cross connections between pores was a 
major structural feature of the porous media (Fatt, 1956). The sphere pack model did 
have interconnected pores, but the shapes of these pores were so complex that no 
analysis of flow through them was possible. If the two models are combined by 
substituting a uniform cylindrical tube for each pore space in the sphere pack model, a 
three dimensional network of tubes is obtained on which, in principle, exact flow 
calculations can be made (Fatt, 1956). 

2.2 Network Model 

Network models represent an interconnected network of capillary tubes within a porous 
media. The lattice of cubic chambers and rectangular tubes of the network represent 
the pore bodies and the pore throats respectively. These models can be used to 
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construct a picture of the microscopic phenomena in a porous medium and serve as a 
tool for analyzing experimental and theoretical concepts. 

The use of network model is not new in the world of modeling of flow through porous 
media. Such models are used in petroleum engineering, chemical engineering and 
physics and increasingly in hydrology (Reeves and Celia, 1996). These models are 
now heavily used as investigative tools to study the nature of fluid flow from the 
perspective of the pore scale. Majority of the studies have been limited to multi-phase 
fluid flow in porous media. In groundwater hydrology, many such applications were 
reported for pressure-saturation relationship for unsaturated fluid flow (Reeves and 
Celia, 1996; Held and Celia, 2001). However, no study so far has been reported for the 
particle transport through porous media under saturated condition using the network 
models. 

Network models of pore structure were first pioneered by Fatt (1956) based on the idea 
that pore space may be represented as an interconnected network of capillary tubes 
whose radii represent the dimension of pores within a porous medium. One of the early 
attempts was the work of Chatzis and Dullien (1977), who applied bond percolation to 
regular two-dimensional and three-dimensional lattices. Several attempts to develop 
porous media models have been reviewed (Payatakes and Dias, 1984; Dullien and 
Batra, 1970). 

To simplify the network operations, it is desirable to replace the three-dimensional 
network by a two-dimensional network. In terms of the porous medium the 
replacement is justified by assuming that a very thin slice of porous medium with 
impermeable planes sealing the two large surfaces has the same properties as a cube of 
the same material. This is equivalent to assuming that the change produced by making 
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the porous medium thin, and thei*eby eliminating a number of cross-connections in the 
third dimension, may be compensated by introducing additional channels within the 
two-dimensional network. It is also assumed that changes in the tube radii distribution, 
network configuration and other variables influence the properties in the same 
qualitative manner in the two dimensional network as in the three dimensional network 
(Fatt, 1956). For single-phase flow, “blind” or “dead end” pores in the porous medium 
are assumed to be absent. There is very little experimental information on the existence 
of blind pores in real porous media. 

2.2. 1 Spatial Distribution of Tubes 

The network model generally chosen to represent the porous media is a regular network 

lattice having nodes in X and Y directions respectively. The nodes of lattice represent 

/ 

the pore bodies; each node is connected with pore throats, which are represented by 
bond of the lattice. The only parameter used to describe the network configuration is 
the number of tubes connected to each tube. This parameter is typically called the 13- 
factor. Many different regular networks can be chosen (Fatt, 1956). Some of these are 
shown in Figure 2.1 and their (3-factors are tabulated in Table 2.1. 

Table 2. 1 3-factors for four networks shown in Figure 2.1. 


Network 

3 Factor 

Single hexagonal network 

4 

Square network 

6 

Double hexagonal network 

7 

Triple hexagonal network 

10 


- 6 - 





Chapter2 


Literature Review 



Square Network 


Hexagonal Network 



Triple Hexagonal Network Double Hexagonal Network 

Figure 2.1 Different types of network used for the description of pore network in 
porous media (Fatt, 1956). 

2.2.2 Pore Size Distribution 



a 


To this end the most realistic models for porous media are those which approximate the 
pore space as a network of pores in which relatively larger pores (pore bodies) are 
connected by relatively narrow pores (pore throats). Existing pore size determination 
techniques (eg. Mercury porosimetry) provide information only about an equivalent 
capillary size, which is heavily weighted towards the smallest pore throat dimension 
(Jenkins and Rao, 1984; Dullian and Dhawan, 1974). Ioannidis and Chatzis (1993) 
presented a methodology for estimating pore throat and pore body size distribution of a 
large class of porous media for use in the network model. For the network model to be 
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geometrically defined, dimensions of all pore bodies and pore throats must be 
determined. It is assumed that the pore throat and pore body characteristic dimensions, 
b, and bb follow Weibull distributions, f, and fb respectively. These are generally 
defined as follows (Ioannis and Chatzis, 1993): 


/,(*/)'= 


(A] 


UJ 

«» J 


F, (£>,), i = b,t 


(2.1) 


m) 


exp 

f b i b i min 

y v 

-exp 

r fb. ~b. . y i 

_ l a , . 

_ l J _ 

l-exp 

f 

V 

max min 

a i J _ 



( 2 . 2 ) 


where, 


b, mi,, = minimum pore throat radius, 
b/ max = maximum pore throat radius, 
b, = random number obtained from the distribution, 
<2i and = parameters of distribution. 


Once the frequency distributions/, and f, are selected, pore sizes b, and bb are assigned 
to the nodes and bonds of the network according to the bond-correlated site percolation 
scheme (Diaz et al., 1987). The pore throat, are not distributed randomly over the 
network, but are correlated to the pore bodies that are randomly distributed. This 
correlation is expressed by the following equation (Chatzis, 1980): 

F,(b,) = F b 2 (b b ). (2.3) 

The physical significance of this correlation is that, large pore throats are associated 
with the large pore bodies. Such correlations have been experimentally observed by 
Wardlaw et. al. (1979). In the construction of the network model the node-to-node 
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distance (/) is adjusted so that the porosity of the network is equal to the porosity of the 
porous media under consideration. 

2.2.3 Simulations Using Network Model 

Network models of pore structure, based on bond-con-elated site percolation concepts, 
were successfully used to predict the drainage capillary pressure curves for mercury-air 
and oil-water displacements, absolute permeabilities and formation resistivity factors of 
sandstones (Ioannidis and Chatzis, 1993). Reeves and Celia, 1996 used pore scale 
network model to determine quantitatively a macroscopic relationship between 
capillary pressure, fluid saturations, and interfacial area between fluid phases. Network 
model have been used to study diffusion and dispersion (Hollewand and Gladden, 
1992; Burganos and Payatakes, 1992; Koplic et al., 1988), pore-scale evaporation 
process (Nowicki et al., 1992), the formation and flow of foam (Laidlaw et al., 1993), 
the interpretation of mercury porosimetry and associated characterization of pore size 
distributions (Ionnidis and Chatzis, 1993; Chatzis and Dullien, 1985), and ganglion 
formation . and mobilizations (Dias and Payatakes, 1986). There have been several 
network model developed to study the constitutive relationships required to 
parameterize multiphase flow equations. The relative permeability-saturation function 
has been studied, primarily by petroleum engineers (Bryant et al., 1993; Blunt and 
King, 1991; Kantzas and Chatzis, 1988). The moisture retention function for both two 
and three phase systems has also received attention (Soli and Celia, 1993; Soil, 1991; 
Ferrand and Celia, 1989). 

2.3 Bacteria Transport Models 

Porous media refer to a fixed bed of granular material containing the pores through 
which a solution may flow. When flowing through porous medium, the particles are 
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brought into contact with the particle retention sites; they are retained there or are 


carried away by the stream. Particles that collect on the porous media form a deposit 


that can alter the fluid flow properties and decrease the permeability. 


The mechanisms by which particles in the suspension are removed within a porous 
media are complex and can be described in three steps: (i) Transport of the particles to 
the collector surface or mass transfer efficiency ( 7 ), (ii) attachment of the particle to the 
collector surface or adhesion efficiency (a), and, (iii) detachment of particles from 
collector surface due to hydrodynamic forces. Transport and adhesion are physical 
processes, which depend on the characteristics of suspended solids (size, shape and 
density), characteristics of collector (surface area, pore size and shape) and 
characteristics of flow (velocity, viscosity and temperature). 


2.3. 1 Models Based on Colloid Filtration Theory 

The fundamental equation describing the particle retention in a filter bed is the 
macroscopic mass conservation equation. The basic filtration equation proposed by 
Iwasaki (1937) for the particle removal from the solution is as follows: 

— = -A 0 C with C = C 0 at x = 0 (influent) (2.4) 

dx 

where, 

Ao = initial filter bed coefficient. 

Solution of Equation 2.4 is, 

C{x) = C 0 exp (~A 0 x) 

The equation is only valid for a clean bed, i.e. at the initial stages of filtration. 
Accordingly, Ao is termed as clean bed filter performance coefficient. Initial attempts 
at relating Ao to particle, media, and flow properties were widely divergent since 


- 10 - 



Chapter2 


Literature Review 


usually only suspension concentration was measured and influent particle size was 
poorly known (Ives, 1970). A mechanistic understanding of clean bed filtration has 
followed advances in the modeling literature (Yao et al 1971; Spielman and FitzPatrick 
1973; Payatakes et al. 1974; Rajagopalan and Tien 1976; Tien and Payatakes 1979). 
These models were primarily based on Happel’s spehere in cell model that substitutes a 
single isolated collector for porous media with a liquid shell around it (Happel 1958). 
Deposition is controlled by two processes: transfer of particles from bulk water to the 
collector surface and adhesion. Colloid filtration theory was used to quantify these two 
processes in terms of the collector mass transfer efficiency X (the probability of 
particle approaching a collector to reach its surface) and the adhesion or collision 
efficiency a (the probability of a particle to attach upon reaching the surface). The 
trajectory of a particle near the collector surface was obtained by considering various 
surface and body forces such as, inertial force, gravitational force, drag forces, DLVO 
forces etc. The mass transfer efficiency or the initial collection efficiency (rj) was 
computed based on the assumption that all the particles that reach the collector surface 
are attached, i.e., a- 1. This initial collection efficiency was approximated well by the 
following semi-empirical equation (Rajagopalan and Tien, 1976): 

77 = 4A] / 3 A " 2/3 +0.72A.X;X 5/8 +2.4xIQ- 3 A s N 1 g 2 H; oa ( 2 - 6 ) 

where, 

A Izll 

Av o Q 

2 2 

P = 0-*o) X 

N r , = Wa,/D t 

N, =aja, 
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N G =2al{p p -p)gl{9fOJ) 

_N nht , = A hx(w) / (9na 2 b U ) 

£o = initial porosity of the column, 

U = velocity of the fluid phase, 

D b = diffusion coefficient of the bacterium; 

pb and p = densities of the bacterium and fluid phase respectively, 

a.b = radius of bacterium, 

a s = radius of the media particle, 

p = dynamic viscosity of the fluid phase and 

Abs(w) = Hamaker constant for Van der Waals interaction between a bacterium (b) and 
solid phase (5) across the medium water (w). 

The first term on the right hand side of Equation 2.6 describes the contribution of 
convection and diffusion to 77 while the other terms account for direct interception and 
deviation from trajectories due to other influences. This model for 77 was later 
modified by Basu (2001) to incorporate the effect of hydrophobic forces. The initial 
filter performance coefficient (/to) was related to the initial collection efficiency {a) as 
follows (Rajagopalan and Tien 1976): 

A -I QzlA rfa (2.7) 

4 a x 

A complete mass balance of particles requires consideration of accumulation of 
deposited particles as well as particles in suspension. Assuming uniform flow velocity, 
negligible dispersion, and small local change in suspended solids concentration, the 
mass balance expression becomes (Herzig et al., 1970), 
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where, 

<7= mass deposited per unit volume of bed, 

C= concentration of particles in suspension. 

By further considering the scouring of the deposited mass as a first order reaction, one 


obtains (Herzig, et. al. 1970): 


0(7 

HT 


r, SC 

-U — sc 7 

dx 


(2.9) 


where, 

s = scouring coefficient. 

As particles accumulate, the clean bed filter coefficient Xq no longer applies and 
X changes with the amount of deposition (a). A phenomenological model that relates 
the filter coefficient during clogging of deposited particle volume is as follows (Ives, 
1970), 


X — A 0 




<7y 

^"vmax 


x’ 


( 2 . 10 ) 


where, 

X = filter bed coefficient at time t, 

£= porosity at time t, 

< 7 V m ax - maximum deposit volume per bed volume, 
cr v = volume of deposit per bed volume. 

P,x’,y’ and z’ are positive empirical parameters evaluated from specific filtration data. 
The Equation 2.10 is based on the three processes: an increase in the filter coefficient as 
particle deposit causes an increase in surface area available for future particle 
collection; a decrease in particle filtration as pore spaces are filled and virgin surface 
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area available for particle collection decreases; and decrease of particle collection as the 
fluid velocity within the pores increases due to constriction caused by deposition. 

Alternatively, the adhesion efficiency ( cc ) can be modified as the deposition progresses 
to reflect the change in the filter coefficient A. The level of cc is controlled by cell-solid 
interactions and by the amount of previously attached particle (bacteria). An attached 
bacterium can reduce the deposition by blocking a part of the collector surface. An 
expression for the influence of cell-solid interactions and cell-cell interactions 
(blocking) on a was given by Rijnaarts et al. (1996) 

a = a 0 (l-B6) ( 2 . 11 ) 

where 

ao- clean bed collision efficiency, 

B = blocking factor 

6 - fraction of surface covered, i.e. the number of the cells. 

Initially, 6=0 and a is solely determined by cell-solid interactions; £%= 1 when these 
interactions do not inhibit the adhesion step. B accounts for the screening of the solid 
surface by attached cells, i.e. it is the ratio of the blocked area to the geometric area of 
the cell. It is related to the maximum surface coverage 6 max for single layer adhesion: 
6 m(a = 1 IB. The B factor may vary between 1.5 to more than 20 and depends on 
geometric or hydrodynamic parameters a h and N pe and particle-particle interactions 
(Dabros, 1993 and Admczyk et al, 1992). 

2.3.2 Large Scale Bacterial Transport Models 

Attempts have been made (Corapcioglu and Haridas, 1985; Harvey and Garabedian, 
1991; Homberger et al., 1992) to model microbial transport through subsurface aquifers 
with varying degree of success. All of these models use advection-dispersion-reaction 
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equation where reaction term was appropriately modified to account for deposition, 
detachment, growth and death of bacteria. Deposition and detachment were mostly 
incorporated in a similar form as the mass balance equations described above 
(Equations 2.8 and 2.9). Growth and death were modeled using Monod kinetics and 
first order endogenous respiration, respectively. 

2.3.3 Bacterial Transport Experiments 

In order to investigate the various mechanisms affecting bacterial attachment, a number 
of laboratory studies, involving passage of bacterial solutions through columns packed 
with sand. Teflon grains, glass beads etc., were conducted (Rijnaarts et al., 1996 a and 
b; Fontes et ah, 1991; Scholl and Harvey, 1992; Martin et ah, 1992; Smith et ah, 1985; 
Gannon et ah, 1991). The parameters studied were the size distribution and surface 
properties of collectors (Sand, Teflon and Glass beads etc.), size and surface properties 
of bacteria and, the pH and ionic strength of water. 

An increase in grain size resulted in a decrease of total surface area and hence the 
bacterial adhesion (Fontes et ah, 1991). Bacteria with negatively charged surfaces were 
adsorbed more when the adsorbent was made positively charged (Scholl and Harvey, 
1992). Increase in the ionic strength of the groundwater increased the attachment by 
double layer compression (Martin et al., 1992; Fontes et ah, 1991). Studies at different 
pHs showed that bacteria retention was higher at lower pH and decreased with 
increasing pH (Goldschmid et ah, 1972; Scholl and Harvey, 1992). Gannon et ah, 1991 
studied the effect of the size and surface properties on their transport. Bacterial 
attachment was statistically related to cell size, with bacteria shorter than 1 pm usually 
showing higher percentages of cells being transported through the soil. Othei important 
parameters governing the attachment and transport of bacteria were reported to be 
surface charge and hydrophobicity (Sharma et ah, 1985; Rijnaarts et ah, 1996). 
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From the literature survey, it is clear that the transport of bacteria through the porous 
media is controlled by many different physical parameters such as, advection, 
dispersion, diffusion, size and heterogeneity of pore spaces etc. and body forces such 
as, gravitation, buoyancy, hydrodynamic drag, Van der Waals interaction, electrostatic 
interaction, hydrophobic forces etc. It is very difficult to perform experiments with 
porous media to study the effect of each of the individual parameters on the bacteria 
transport process. In addition, it is difficult to replicate these experiments exactly as the 
pore space geometry, tortuosity and flow condition will vary each time the porous 
media is packed anew. It is also difficult to alter one of these interactions 
independently in the experiments to study the effect, e.g. altering surface characteristics 
for electrostatic interaction will also alter the hydrophobic interaction. If a computer 
model can simulate the characteristics and physics of bacterial transport through porous 
media, some of these experiments can be replaced by simulation or the model can be 
used to design experiments in a logical way. With this objective in mind, the goal of 
the present work was set as follows: 

• Construct a porous media network model for transport of colloid under 
saturated flow condition. 

• A numerical simulation study for bacteria transport in the porous media using 
the network model. 

• Numerical simulations to study the effect of various properties of porous media 
and different hydraulic conditions on the breakthrough of bacterial transport 
through a column. 
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4. 1 Network Model 

4. 1. 1 Spatial Network 

The square network was chosen for the present study because it is relatively easier to 
implement flow and transport models on a square network. However same framework 
can be used to model other type of networks. The nodes of the network and the tubes 
connecting them were given the position numbers separately starting from the lower 
left-hand comer (Figure 4.1). The flow in the network is from bottom to top. Size of 
the network was m x n where, m is the number of nodes in the horizontal direction 
(cross section of the column) and n is the number of nodes in the vertical direction 
(length of the column). 


4. 1.2 Pore Throat (Tubes) Radius Distribution 

Radius of each pore throat in the network was determined using truncated Weibull 
distribution function restricted in (b, min , b t max ). Cumulative density function of 
truncated Weibull distribution function is given as follows (Ioannidis and Chtazis, 


1993): 
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b t min = minimum pore throat radius, 
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Figure 4.1 Number sequence scheme for the nodes and the tubes of the square 
network. 

h n = number of horizontal nodes 
v„ = number of vertical nodes 

i = number of horizontal layers in the network column 
k= ih„ 
r= 2h n -\ 

T\=ir+(h„-l) 

Total number of nodes = h n v n 

Total number of tubes = h„ ( v„-l)+v„ (h„-\) 
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£Vmax — maximum pore throat radius, 

b t — random number obtained from the distribution and 

a, and # = parameters of Weibull distribution. 

A random number generator (Hull and Dobell, 1962) was used to obtain the random 
numbers from a uniform distribution U(0,1). The random numbers obtained from the 
uniform distribution were transformed to the truncated Weibull distribution using the 
following transformation equation: 

b , = b , min + a, [- In (R v (l - A) + A)]/a 


where, 

Ru = random number obtained from the uniform distribution in (0,1) and 


A = exp 


b -b . 

t max min 


V 1 


Thus, a set of independent random numbers having a distribution given in Equation 4.1 
was obtained. A total of N t (number of tubes in the network) random numbers were 
generated and was assigned sequentially to the tube numbers in the network (Figure 
4.1). This way, the radii of the tubes were spatially randomly distributed in the 
network. The ranges of different parameters of distribution used in this study are 
shown in Table 4.1. 

Table 4.1 Input to get Radius of the Tubes 


Parameters 

Values 

References 


1.0 

Ioannidis, 1993 

A 

2.0 

Ioannidis, 1993 

bt min 

25 pm 

Wan and Wilson, 1994 

bt max 

50 pm 

Wan and Wilson, 1994 
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4. 1.3 Pore Throat (Tubes) Length 

Lengths of the tubes were constant in the regular network. Depending on the soil type, 
either the length of the pore throat or the porosity of the network was supplied. If the 
porosity of the network was supplied as input then the pore throat length was calculated 
and if the length was supplied then the porosity of the network was calculated. The 
calculation was based on the aerial porosity as a two-dimensional network was assumed 
to be similar to an interface, a limiting case of a three dimensional network when one of 
the dimension approaches zero (Hassanizadeh and Gray, 1979). The mathematical 
relation between length of the pore throat (tubes) and the porosity of the network was 
as follows: 


£ ( k n ~ 0 ( V „ - 0 

where, 

l = length of one tube, 

R, = radius of zth tube, 

£ = porosity of the network, 
h n = number of horizontal nodes, 
v„ = number of vertical nodes, 

N t = total number of tubes in the network. 


(4.3) 


4.2 Flow Through the Network 

4.2. 1 Pressure Head at Nodes 

The head loss was calculated at each tube of the square network model using Darcy- 
Weisbach equation. 

h (4.4) 

" ' 4 
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where, 

hki = head loss in the tube connecting nodes k and i, 
vu = velocity in the tube connecting nodes k and i, 
f= friction factor, 

Rki= radius of the tube connecting nodes k and i. 

The above equation for a single tube connected between the nodes k and i can be 
written as: 

K K„ .*/**(&, Ha, p (4.5) 

where, 

Kki = head loss coefficient for the tube connecting nodes k and i, = — ~ — — for n ki - 2 

4 x 2 gR 5 kl 

Hi = head at the ith node, 

Qki= discharge in the tube connecting the nodes k and i, 

nki = exponent in the head loss equation for the tube connecting nodes k and i, 


Since the head loss is positive in the direction of flow, sign of the head loss is same as 
sign of the flow. So the Equation 4.5 can be inverted as follows: 


Q ki =sign{H k -H,] 


Hk-H, 


Jt*hi 


K 


(4.6) 


u j 


The continuity equations for the node i in the network can be written as follows: 

ntf ' 

Z^Qk, = U, for f= 1,2, 3 N, (4-7) 

*=t 

where, 

Ui = consumptive use at node i, 
mi = number of connections to node i. 
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The energy and continuity equations can now be combined as follows taking U ; — 0 for 


each node i : 


- H, 


k = 1 


H k -H, 
K ki ' 


\ Ut hi 


= 0 


for i = 1, 2, 3 N, 


(4.8) 


There is one such equation for each node and the variables are heads (Hi s) at each 
node. So, this gives a system of non-linear simultaneous equation of size h n v„ -2 h„, 
where h n v n is the total number of nodes in the network. This system of non linear 
equations were solved using Newton-Raphson method. Only the heads at the inlet 
boundary and the outlet boundary were supplied. This is equivalent to supplying the 
difference in head between inlet and outlet of a column or the pressure drop in a 
column. / 

4.2.2 Discharge in the Tubes 

The net discharge <2w in each tube connecting nodes k and node i was calculated from 
the Equation 4.6. 

4.2.3 Velocity of Fluid in Tubes 

The net velocity of fluid in each tube was calculated from the net discharge in the tube 
and the cross-sectional area of the tube. 


where, 

Vki = Flow in the tubes connecting the nodes k and i. 
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4.3 Bacteria Transport Model 

4.3.1 Inflow Concentration 

After determining the velocity of the flow in all the tubes of the network the influent 
concentration Co = 0.005 kg/m 3 was added to the fluid going into the network. In all 
the tubes placed at the bottom layer of the network having the position number h n to 
2/in-l were given an initial concentration Co = 0.005 kg/m 3 and the rest of the tubes 
were having Co = 0 at time t = 0. It may be noted that in this study, SI system was used 
for unit of all the parameters. If the influent bacteria concentration is defined 
alternatively as the number of cells per 100 mL, the values and units of all other 
parameters need to be changed accordingly. 

4.3.2 Effluent Concentration and Bacterial Deposition 

The tubes in the network are very small in length. For small time step (At), A may be 
considered independent of x and t in each tube during every time step. This enables the 
Iwasaki equation (Equations 2.4 and 2.5) to be applied to individual tubes for the 
calculation of effluent concentration. In the Equation 2.5, the Ao was substituted by A 
and it was modified after every time step using Equation 2.5. The influent 
concentrations (Co,) were the concentrations at the influent nodes of the tubes and 
changed at every time step depending on the concentrations at various nodes. 

For the calculation of < 7 , detachment or scouring can be neglected (5 = 0) as bacterial 
deposition is often considered to be irreversible (Rijnaarts et al. 1995, Sjollema et al. 
1990). The deposition ( 0 ) can be obtained by combining the Iwasaki Equation 2.4, and 
Equation 2.8 and by substituting A for Ao to take into account the time dependence. It is 
possible to obtain a closed form solution of Equation 2.8 for a small time step in a small 
tube (constant A), as follows (McDowell-Boyer, 1986): 
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Due to deposition, the porosity of the network changes. Hence the velocity was also 
modified at each step using the relation: 

v = v 0 (l-<7) (4.11) 

The maximum deposition was taken as 10% in the present study. 

The effluent concentration was calculated at the end of each time step (At) from the 
individual tube using the Equation 2.5 with the updated value of A and C 0 is taken as 
the inflow concentration in the tube at that time. The final effluent concentration from 
the column of the network was net concentration of out flow from the vertical tubes at 
the upper end. 


4.3.3 Filter Coefficient 

The initial filter coefficient ( Ao ) was calculated for each tube using the Equation 2.7. 
The values of the filter coefficients (2) were modified using Equation 2.10 after every 
time step (At) for each tube. 


4.4 Numerical Scheme of Solution for Breakthrough Curve of 
Bacterial Transport 

Flow chart of the whole network model is shown in Figure 4.2. A detailed outline of 
the time stepping scheme involved in the calculation of effluent concentration (C) and 

deposition (o) is given below: 
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Step 1: The minimum time taken by the fluid coming from the bottom layer, to reach 
each of the nodes of the network was calculated. Hence each node was associated with 
the minimum time T m j n . 

Step 2: To determine the concentration profile the time step At was chosen. 

Step 3: Concentration of the particle was calculated at each node after the time step At. 
If T m in at any node is greater than At, then concentration at that node was taken as zero. 
Step 4: After each time step At, the total deposition in the tube was calculated, the 
velocity of the fluid in the tubes was modified due to the bacterial deposition and the 
new filter coefficient (A) was calculated. 

Step 5: At each node the concentration was determined and effluent concentration was 
calculated. 
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Figure 4.2 Bow chart representing the steps involved in bacterial transport model 
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Results and Discussion 

5.1 Generation of Porous Media Network 

5. 1. 1 Effect of a t on the Radius Distribution 

The effect of C(\ on the radius distribution is shown in the Figures 5.1a-b. For a 
constant value of y 3, = 2.0, the radii distributions were compared for a, = 0. 1 and a, = 1. 
The shape of the distribution curve was observed to remain same in both the cases but it 
was shifted towards higher radius values. 

5. 1.2 Effect of ft on the Radius Distribution 

/■ 

The effect of ft on the radius distribution is shown in the Figures 5.2a-b. For a constant 
value of ^ =1.0, the radii distributions were compared for ft = 1.0 and ft = 4.0. The 
values of ft are related to the shape of the distribution curve. For the higher values of 
ft (= 4.0) the distribution was more symmetric and narrower compared to that with the 
lower value of ft (= 1 .0). 

5. 1.3 Effect of b t max and b t m m on the Radius Distribution 

The parameters b, max and b, m j n determine the mean pore size as well as the range of the 
distribution and are directly related to the soil type. For example, sandy soil is likely to 
have bigger pore sizes compared to those of clay soil. Three different ranges of pore 
size distribution with same <3 and ft values are shown in Figures 5.1b and 5.3a-b. 
These distributions were generated by changing the values of b, max and b, m >n- 
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5. 1.4 Validation of Porous Media Network through Pressure Saturation 
Curve 

Generation of pressure saturation curve does not require solving the hydraulic model or 
the tiansport model. It can be directly obtained from the network properties (Fatt, 
1956). For 

the present study, the value of cx t and j3, were taken as 1.0 and 2.0, respectively. The 
relation between capillary pressure and wetting phase saturation is shown in Figure 5.4. 
The network data gives the stepped curve because the radius distribution is not 
continuous. In real porous media with continuous pore throat radius distribution, the 
capillary pressure saturation curve is always smooth. In Figure 5.4, the dashed line 
shows the result of smoothing the network data. The relation obtained from the 
network model resembles a typical capillary pressure curve with the exception that the 
curve derived from the network is asymptotic to zero wetting phase saturation, whereas 
curve from real porous media is asymptotic to 20 to 30 percent wetting phase saturation 
(residual saturation). This difference is due to the dead end pores in the real porous 
media that can trap fluid during drainage. However, the curve matches with other 
reported results from network model (Fatt, 1956). 

5.2 Hydraulic Model 

The hydraulic model was tested with small network sizes of 3x3 and 3x4. The values 
obtained from the program were compared to the manually calculated values. It was 
tested for both constant tube radii and radii distribution obtained from Weibull 
distribution. The simulated heads at the nodes and the corresponding discharges in the 
pipes matched the computed values exactly. Figure 5.5a and b shows a typical result. 
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5.3 Bacterial Transport Model 

Following the confirmation of the porous media network model and hydraulic model, 
the bacterial transport model was incorporated in the network. The following sections 
describe the verification of the bacterial transport model for consistency. 

Table 5.1 Input supplied 


Parameters 

Notation 

Values 

References 

Particle diameter 

dp 

lx 10‘ 6 m 

Basu, 2001 

Grain diameter 

dm 

0.54x1 O' 3 m 

Basu, 2001 

Relative particle density 

Pp 

1.05 

Basu, 2001 

Fluid density 

Pf 

1.0 x 10 3 kg/m 3 

Basu, 2001 

Friction factor 

f ' 

0.05 

Subramanya 2000 

Diffusivity 

D P 

4.3x 10‘ l3 m 2 /sec 

McDowell-Boyer, 

1986 

Initial Porosity 

£o 

0.4 

Basu, 2001 

Initial bacteria 

concentration 

Co 

0.005 kg/m 3 

Basu, 2001 

Fluid dynamic viscosity 

M 

0.00098 N-s/m 2 

McDowell-Boyer, 

1986 

Empirical parameters 

x, y, z, 0 

1.0 



5.3. 1 Effect of Time Step on the Breakthrough Curve 

The choice of time step plays a critical role in the transport model. If the time step 
chosen is larger than the smallest flow through time in one of the tube of the network, 
deposition of particles will be time-lumped in more than one tube. The effect of time 
step on the breakthrough curve obtained from the network is shown in Figure 5.6. For 
this simulation, a network size of 10x10 and conservative tracer (or = 0 in Equation 2.7) 
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was used. The breakthrough curves became invariant of the time step beyond 2.5xl0’ 6 
sec. Based on this result, the time step was taken as 2.5x1 0" 5 sec for all the simulations. 

5.3.2 Effect of Cross Section of the Column 

The breakthrough curve is expected to be invariant of the column cross sectional area. 
This was checked by keeping the number of vertical nodes fixed at 10 and expanding 
the network horizontally by increasing the number of nodes from 10 to 50. The 
simulations were once again conducted for a conservative tracer with zero retention. 
Figure 5.7a shows the results of these simulations. As expected, the breakthrough 
curves remain invariant of the column cross section. The network model imparts 
dispersion in the breakthrough curve (Figure 5.7a) evident from the gradual variation of 
the effluent concentration before it reaches the influent concentration. This is typically 
observed in the experiments (Fontes et al., 1991; Harvey and Garabedian 1991). 
Simple application of the Iwasaki Equation (2.4) to a single tubular column reactor will 
give no dispersion. This indicates that a network model is capable of producing results 
that is characteristically similar to the experimental results. 

5.3.3 Effect of Length of the Column 

Increasing the length of the column is likely to delay the initial breakthrough as well as 
impart more dispersion (Fontes et al., 1991; Harvey and Garabedian 1991). In the 
simulations, the number of horizontal nodes was kept fixed as 10 and the network was 
expanded vertically by increasing the number of nodes from 10 to 50. Once again, the 
simulations were conducted on a conservative tracer. The effects of increasing the 
number of nodes in the vertical direction on the breakthrough curve are shown in the 
Figure 5.7b. As expected, the time for initial breakthrough as well as the dispersion 
increased with increasing length of the column. 
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5.4 Comparison of Simulated Results with Experimental 

Results 

A set of experiments was reported (Rijnaarls et al. 19%) for bacterial transport through 

same media but different lengths of column. The model simulation was conducted to 

see how a 2-d network model simulates the results obtained from 3-d column 

experiments. The pore size distribution as well as the adhesion efficiency ( n ) was 

unknown. The pore size distribution was obtained by taking o; =1.0 and $ = 2.0. The 

adhesion efficiency was assumed as unity. Simulation results (Figure 5.8a) were 

compared to the experimental results (Figure 5.8b). The time scale of breakthrough 

was much smaller in the 2-d network than the experiments as the pore volume of the 

network is much less and an added dimension in the columns gives the pores more 

* 

connectivity and tortuosity. However, similarity of the effect of column length on the 
breakthrough levels of C/Co is evident. The point at which C/Co reached a plateau 
decreased with increasing column length similarly in simulations and experiments. 

5.5 Effect of Different Physical Parameters 

5.5. 1 Effect of Porosity 

In order to study the effect of porosity, a network of dimension 10x100 was taken and 
the porosity was changed keeping all other parameters of the model constant. The 
comparison of (he breakthrough curves for three different porosities 0.3, 0.4 and 0.5, is 
shown in Figure 5.9. It shows that for the same soil type and bacteria, a well- 
consolidated bed of soil will transport higher amount of bacteria due to lower pore 
space available'Tor deposition. 

5.5.2 Effect of Flow Rate 

For the same column, lower flow rate yields higher hydraulic retention time. This is 
equivalent to increasing the length of the column for the same flow rate. So, the effect 
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of decreasing the flow rate is same as increasing the length of the column and is 
demonstrated in Figure 5.10 for a 10 x 100 network. 

5.5.3 Effect of (Tv max 

Another parameter that is very hard to determine experimentally is the maximum 
specific deposit (<j v max) beyond which the no more bacterial deposition is possible. As 
a result, it is impossible to set a value of < 7 V max for simulation. This set of simulation 
was conducted to see the effect of the value of maximum specific deposit on the final 
breakthrough curve. It is easily seen that the variation, although present, is not very 
significant (Figure 5.11). A decrease of one order of magnitude to the value of <r v max 
resulted in less than 10% difference in the breakthrough curve. Although, it is expected 
that, at the lower values of cr v max , the result will be more sensitive since, X approach 
infinity as <r vmax approach zero Equation 2.10. 

5.5.4 Effect of Adhesion Efficiency (a) and Mass Transfer Efficiency ( rj) 

At present, no model or description exists for the adhesion efficiency (a). It is often 
chosen arbitrarily. Rijnaarts et al. (1996) proposed various ways to compute the 
adhesion efficiency from the breakthrough curves of column experiments but the 
results varied depending on the method used. As a result, it is impossible to set this 
parameter for any simulation. A sensitivity study was conducted to observe the effect 
of cron the breakthrough curve (Figure 5.12). For the full range of a (0 to 1), the 
breakthrough curves vary significantly indicating the adhesion efficiency to be an 
important parameter. Since, the mass transfer efficiency and the adhesion efficiency 
are multiplied to each other in the expression of the filter coefficient (Equation 2.7), the 
breakthrough curves are equally sensitive to the mass transfer efficiency. 
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5.6 Application Perspective 

Application of the 2-dimensional network model is in designing column experiments 
for the colloid transport through a packed bed or porous media. The model is capable 
of predicting trends as well as variations of the breakthrough curves under different 
parameters and conditions. This will help in deciding the key control variables and 
parameters for the design of colloid transport experiments in porous media. Thus, it 
will reduce the number of experiments to be conducted. Expanding the 2-d model to a 
3-d network is possible but the computational requirement may outweigh the advantage 
of conducting less number of experiments. The change from a 2-d to a 3-d network is 
unlikely to provide any added insight into the mechanism of problem. Lastly, 
calibration of the network model parameters with experimental data may enable the 
model to be used as a predictive tool. 


- 33 - 




Radius, p m 


Figure 5.1a Effect of distribution parameters on the radius distribution. b lmux = 50.0 
pm, btm-m = 25.0 pm, a,= 0. 1 , J3,= 2.0. 
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Figure 5.1b Effect of distribution parameters on the radius distribution. b tm(tx = 50.0 
pm, btmm = 25.0 pm, a,- 1.0, A = 2.0. 
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Figure 5.2a Effect of distribution parameters on the radius distribution. = 25.0 
P-ifi, b lm(IX = 50.0 ptm, a,= 1.0, #=1.0. 
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Figure 5.2b Effect of distribution parameters on the radius distribution. b, min = 25.0 
|im, b t max= 50.0 |im, a,=\ .0, ^,=4.0 
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Figure 5.4 Nondimensional capillary pressure characteristics of square network of 
size 10x10 using radius distribution of Figure 5.1b. 
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Figure 5.5a Computed values of pressure head for 3x3 Network. Black numbers 
show the supplied values for the head at the end nodes. Green numbers 
are computed head at the nodes for the uniform radius network. Blue 
numbers represent computed head at the nodes for the network having 
radius distribution obtained from the weibull distribution. 
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Figure 5.5b Computed values of pressure head for 3x4 Network. Black numbers 
show the supplied values for the head at the end nodes. Green numbers 
are computed head at the nodes for the uniform radius network. Blue 
numbers represent computed head at the nodes for the network having 
radius distribution obtained from the weibull distribution. 
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Figure 5.6 Effect of time step on the breakthrough curve for bacterial transport. 

Network dimension 10x10, Head =1 m, Porosity = 0.4, a =0. Legend 
shows time step in seconds. 
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Figure 5.7a Effect of cross section of the network on the breakthrough curve for 
conservative tracer. Head =1 m, Porosity = 0.4, <2=0, Time step = 
2.5E-5. 



0.01 0.02 0.03 0.04 


Time, Seconds 

Figure 5.7b Effect of vertical length of the network on the breakthrough curve for 
conservative tracer. Head =1 m, porosity = 0.4, <2=0, Time step - 
2.5E-5. 
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Figure 5.8a Effect of vertical length of the network on the breakthrough curve for 
bacterial transport. Head=l m, porosity 0.4, a- 1, Time step 2.5E-5. 



Figure 5.8b Effect of column length on the breakthrough curves for bacterial 
transport. (Rijnaarts et al., 1996.) 
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Figure 5.9 Effect of porosity on the breakthrough curve for bacterial transport. 

Head 1 m, Network dimension 10x100, a = 1.0, Time step 2.5E-5. 
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Figure 5.10 Effect of flow rate on the breakthrough curve for bacterial transport. 

Network dimension 10x100, a= 1.0, Time step 2.5E-5, Porosity = 0.4. 
Legen shows flow in the network. 
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Figure 5.11 Effect of maximum deposition on the breakthrough curve for bacterial 
transport. Head = lm, Network dimension 10x100, a- 1.0, Time step 
2.5E-5, Porosity = 0.4. 
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Figure 5. 12 Effect of oron the breakthrough curve for bacterial transport. Head =1 
m. Network dimension 10x100, Porosity = 0.4, Time step 2.5E-5. 
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Chapter 6 

Summary and Conclusion 

The present research was undertaken in two steps. In the first step the two-dimensional 
network model was developed to represent the porous media. In the second step the 
bacterial transport study was done. The Weibull distribution function was used to 
determine the radius of the tubes of the network representing porous media. The 
network was cestea for successful generation of pressure saturation curve. After 
validating the hydraulic model with the hand-computed values in simple network, the 
model was used to generate the breakthrough curves for bacterial transport. The 
following observations can be made from the study: 

• The distribution parameters, a, and /i, play important roles in determining the 
shape and spread of the pore size distribution. The parameters b x m in and b x max 
set the range of distribution and their values should be governed by the 
knowledge of the properties of porous media. 

• Choice of time step is an important factor for the consistency of the model 
results. 

• Application of Iwasaki equation over the individual tubes in a network with the 
consideration that the filter coefficient is constant for a small tube in a small 
time step, is capable of producing the expected dispersion trends for a 
conservative tracer. 

• The effect of cross section and the vertical length of the column for a 
conservative tracer produced the expected trends. The breakthrough curve 
remained invariant with the changing cross section of column. Increasing the 
column length increased the dispersion. 
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• For a non-conservative particle (bacteria) transport, the model results simulated 
the trends of the experimental results. 

• In the sensitivity studies, adhesion efficiency (a) was found to be a very 
important parameter for bacterial transport. 
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Chagter 7 

Scope of Future Work 

The Following Areas are recommended for further investigations: 

• A Happel’s sphere in cell model for the mass transfer efficiency (77) may be 
incorporated in the network model. This will enable the model to study the 
effect of different surface forces arising out of various particle and collector 
surface characteristics. 

• A deterministic model for the collection efficiency (a) is needed. 
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